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DIFFERENTIAL EQUATIONS WITH ALMOST 
PERIODIC COEFFICIENTS 
By Thomas J. Coakley 
Ames Research Center 


SUMMARY 


A new method for analyzing linear differential equations with almost 
periodic coefficients is described and discussed. The method is essentially a 
generalization of the classical Floquet theory developed for linear periodic 
systems. By use of the method, it is shown that an almost periodic system 
with rapidly varying coefficients can always be reduced or transformed to a 
corresponding system with either constant or slowly varying coefficients. An 
example is discussed which illustrates the method and indicates how the 
reduced system is frequently more easily analyzed than the original system. 


INTRODUCTION 


Ordinary differential equations constitute one of the principal tools for 
mathematically modeling physical systems. Within this general class, linear 
equations draw special attention because they are tractable and can be used 
frequently to approximate nonlinear systems. Most prominent within the class 
of linear differential equations are those with constant coefficients. These 
types of equations are particularly useful because their theory is quite 
simple and well developed. Unfortunately, the theory of linear differential 
equations with variable coefficients is not so well developed. The present 
investigation is a contribution to this theory and, in particular, is con- 
cerned with linear homogeneous equations with periodic and almost periodic 
coefficients. These equations have important applications in the dynamic sta- 
bility of vehicles, structures and fluids, the theory of orbits, and the 
quantum theory of radiation. 

Compared to the theory of linear differential equations with almost 
periodic coefficients, the theory of linear differential equations with peri- 
odic coefficients is more complete. For these systems, the classical Floquet 
theory (refs. 1 and 2) plays a central role. In this theory, it is shown that 
a linear system with periodic coefficients can always be transformed into a 
corresponding system with constant coefficients. Although Floquet theory is 
useful qualitatively, it is of limited practical application because no algo- 
rithm is provided for the actual construction of the transformation. For this 
purpose, additional and more specialized techniques must be used. 


Ill 


For linear systems with almost periodic coefficients, unfortunately, 
there is no general result such as Floquet theory to relate behavior and sta- 
bility properties to systems with constant coefficients. As a consequence, 
methods for the analysis of these systems are much less developed and more 
highly specialized than those for linear periodic systems. 

In this report, a new method is described which, in a certain sense, 
generalizes Floquet theory to linear equations with almost periodic coeffi- 
cients. Our main result may be stated as follows: Under very broad condi- 

tions a bounded linear transformation can be found by means of which a linear 
system with rapidly varying, almost periodic coefficients may be transformed 
into a corresponding linear system with more slowly varying coefficients. 
Moreover, a useful algorithm is provided for constructing the transformation 
and the transformed system. The utility of the method rests on the fact that 
the transformed system is frequently more easily analyzed than the original 
one . 


ANALYSIS 


The method to be discussed is closely related to the so-called averaging 
methods previously developed for certain classes of periodic and almost per- 
iodic systems (refs. 2-6). However, it is more general and rigorous. This 
permits, on the one hand, application to wider classes of systems and, on the 
other, a more detailed discussion of existence and convergence. 


Basic System 

The basic system to be investigated is 

x = A(t)x (1) 

where x(t) is an n -dimensional vector function, A(t) is a bounded matrix func- 
tion, and x = dx/dt. For the current investigation, it is assumed that A(t) 
can be represented as a generalized Fourier series, or almost periodic func- 
tion, of the form 

A(t) =E(An cos w n t + B sin u) n t) = [a.. (t)] (2) 

n - L J 

where the frequencies, tu n , are arbitrary. The series is assumed to be abso- 
lutely convergent so that the norm of A is bounded, that is, 

II All =E^I A n!l + Kll) < » 

n 

where the matrix norms ll A nll and || B n || are 

H A nll = m f x ? I a nij I etc ' 
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Classical Floquet Theory 


If A(t) is periodic with period T = 2ir/to 0 so that w n = nm 0 in 
equation (2), the classical Floquet theory (refs. 1 and 2) states that a per- 
iodic transformation exists by means of which equation (1) can be reduced to a 
system with constant coefficients. This theory is outlined briefly here. 

Let X(t) be a fundamental matrix of equation (1) so that 

X = A (t) X , X(0) = I (3) 

where I is the identity matrix. It is easily shown that X(t)C, where C 
is an arbitrary nonsingular constant matrix and X(t + T) are also fundamental 
matrices. The latter assertion follows from the periodicity of A(t), that is, 

dT X(t + T) = dcA ' Tj X(t + T) = A(t + T)X(t + T ) = A (t)X(t + T) (4) 
By equating X(t + T) and X(t)C, it follows that 

Now let the matrices B (constant) and P(t) 

e BT = C = X (T) 1 

P(t) = X (t)e -Bt J 

The matrix P(t) is periodic (and nonsingular) since 

P(t + T) = X(t + T)e“ B(t + ^ = X(t)e BT e“ B(t + T) = P(t) (6) 

By making the transformation 

x = P(t)y (7) 

and using equations (3) and (5), it follows that the equation for y is 

y = By (8) 

In carrying out this transformation, several procedures are available. 

If the matrix A(t) is arbitrary, then the reduction can generally only be 
made numerically. This can be done by first solving equation (3) over one per- 
iod, that is, 0 < t < T. The matrix B can then be obtained (often in diago- 
nal form) from equations (5) by solving the matrix eigenvalue problem for 
C = X(T) (refs. 7 and 8). If A(t) contains a small parameter, for example, 

A = A q + eA 1 (t) , where A c is a constant matrix and c is the small 
parameter, then perturbation methods can be used (refs. 2, 4-6). 


X (T) = X(0)C = C. 

(periodic) be defined by 

(5) 
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Averaging and Integrating Operators 


We propose to show that equation (1) can be transformed into a similar 
system but with only low-frequency components. For this purpose, we intro- 
duce two linear functional operators, Mp and Np, called averaging and inte- 
grating operators, respectively, defined as follows: 
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( 9 ) 


The number Up is called the cutoff frequency and is assumed to be 
greater than zero. It is easily verified that the operations M A and N A 
have the properties: ^ — , P P 


M p A 


- L f A n cos u) n t + B n sin to n t) 


°“ a) n <w p 


-L 


) ( 10 ) 


-1, 


p A Z i w n (Aji sm ai n t - B cos ou-t) 

a)- £00 11 11 

P n 


J 


h- V ' A ' V til) 

From these relations, it is seen that MpA, which we choose to call the 
generalized mean or average value of A, contains no frequencies greater than 
u) . Conversely, NpA, which will be called the generalized integral of A, 
only contains frequencies equal to or greater than o) p . It is easily shown 
that MpA and N A are bounded, that is. 


IlMpAll £ ||a|| , 1 1 Np a| ) £ « p “ 1 Ha|| 


( 12 ) 


The operators defined above are closely related to those used by Krylov- 
Bogoliubov and Mitropolski (refs. 3 and 4). In fact, the operators they used 
may be obtained from equations (9) and (10) by letting the cutoff frequency go 
to zero, that is, (Dp -»■ 0. In this case, the function MpA is a constant 

equal to the mean value of A. However, the corresponding integrating opera- 
tion, NpA, is not bounded, and additional and rather severe restrictions must 

be placed on the class of almost periodic functions for this operation to be 
well defined (and bounded) (refs. 5 and 6). 
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Transformation of the Basic System 


The averaging and integrating operators defined above will now be used to 
transform equation (1). For this purpose, we assume that the vector, x, in 
equation (1) is transformed into a new vector, y, by a matrix I + V(t), that 
is, 

x = [I + V(t)]y (13) 

where I is the identity matrix and V(t) is to be determined. Assume fur- 
ther that the transformed differential equation for y is of the form 

y = B(t)y (14) 

where B(t), like V(t), is to be determined. Substituting equations (13) and 
(14) into equation (1) and rearranging, we obtain the following equation for 
V : 

V = U - B (15) 


where 


U = A + AV - VB (16) 

Suppose for the moment that the matrix function U(t) in equation (15) is 
known and of the same form as A(t) (i.e., with an arbitrary frequency spec- 
trum). We desire to choose B(t) in such a way that solutions to equation (15) 
for V are bounded for all t regardless of the form, or frequency spectrum, 
of U(t). It is readily verified that this can be done by taking 

B = MpU (17) 

The solution for V then becomes 

V = N U + V (18) 

p o 

where V Q is an arbitrary constant matrix that will be taken as zero in the 
current development. By means of equation (12), one may verify that V is 
bounded if U is and if o)p > 0. 

In order that B and V, as defined above, yield the transformation given 
by equations (13) and (14), U(t) must satisfy equation (16). By setting 
B = nLu and V = N U, this equation is reduced to a single equation for U, 
that is, ^ 


U = A + A (NpU ) - (N p U)(MpU) 


(19) 


The function U defined by this equation will be called the generating 
function of the transformation. It would be very useful to know that this 
function exists for arbitrarily small cutoff frequencies including the value 
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w p = 0. In this case, the reduced system, equation (14), would have constant 

coefficients, and one would have the direct generalization of Floquet theory 
to almost periodic systems. Unfortunately, because of the unboundedness of 
Np as cop -*■ 0, this is not generally possible. However, solutions to equation 
(19) can be shown to exist for sufficiently large values of the cutoff fre- 
quency, u)p, or for sufficiently small matrices A(t). The explicit condition 
on A and Up which guarantees the existence of solutions to equation (19) is 

|| A|| < (3 - 2/2 )oj p (20) 

If this condition is satisfied, it can also be shown that solutions to 
equation (19) can be constructed by successive approximations, starting with 
the initial approximation U = Uj = A. Further, the matrix NpU is bounded 
by unity, that is, j| NpU|| < 1, so that the transformation matrix I + V of 

equation (13) is bounded and nonsingular. It then follows that the ultimate 
stability and behavior of solutions to equation (1) are controlled by the solu 
tions to equation (14) . These results can be proven by means of the contrac- 
tion mapping principle, or fixed point theorem in Banach space (refs. 5 and 9) 
However, the details of the proof are beyond the scope of the present 
discussion and will be omitted. 

In addition to solving equation (19) by successive approximations, one 
can also use series expansions in terms of a parameter, say, e. In this event 
one writes 


A (t) = eA(t, e) = e n A n (t) 


"\ 


n=l 


( 21 ) 


U(t) = eU(t, e) 


" I>n 


(t) 


n= 1 


where the series for A(t) is assumed to converge in some neighborhood, 

0 < £ < E and || A (t , e)|| < 1. The series for U(t) can be constructed by 

substituting equations (21) into (19) and equating coefficients of e n . It 
can be shown that this series converges if e < e Q < kcop, where k is a 
constant of the same order of magnitude as the parenthetical constant 
appearing in equation (20), that is, 3 - 2/2". 


Some Properties of the Reduced System 

In certain cases, the reduced system, equation (14), may have properties 
that considerably simplify its analysis. This is certainly true when the sys- 
tem matrix, B = MpU, is constant. An important example is the case when A(t) 

is periodic and expressible in the form 

oo 

A(t) = L (A n cos n<jo 0 t + B n sin noo 0 t) (22) 

n=o 
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By choosing the cutoff frequency equal to the base frequency, that is, w = w 0 , 
and assuming the existence criterion, equation (20), to be satisfied, one may 
show that MpA and thus B = MpU are constant matrices. Although this is 

simply a restatement of Floquet theory, there is the additional feature that 
solutions can be constructed by the rather straightforward and convergent pro- 
cesses of successive approximations or series expansions. In this aspect the 
method is similar to methods developed earlier by Cesari (ref. 2), Hale 
(ref. 5), and Golomb (ref. 6). However, in those methods the matrix B (or, 
equivalently, its eigenvalues) is defined implicitly by a nonlinear matrix 
equation that must be solved. In our method, it is defined explicitly, and 
this feature generally simplifies theoretical and numerical analyses. 

In the more general situation in which A(t) is almost periodic, that is, 
where some or all of the frequencies a) n are incommensurable, there is no 
assurance that the matrix B will be constant. In fact, B will almost 
always be time dependent. However, it frequently can be made constant to an 

acceptable order of approximation. An important example is when A(t) can be 

expressed in the series form of equations (21) and each term, ki (t) , A 2 (t), 
etc., is a trigonometric polynomial containing only a finite number of frequen- 
cies. In this case, it is possible to insure that the first N terms in the 
series for B are constant by choosing the cutoff frequency, ajp, sufficiently 
small. When only the first term is to be made constant, that is, N = 1, then 

00 p is chosen equal to the smallest frequency in A* (t) , say, co i . The first 

approximation to B = MpU is then eMpAi = const, and the sufficient condi- 
tion on the size of e for the series to converge is estimated by the condi- 
tion e < ka)i . One may show that as the level of approximation for which the 
matrix B is to be made constant, that is, N, is increased, the required cut- 
off frequency (and thus the sufficient interval of convergence in e) 
decreases so that cop and e 0 in the limit as N The series expansion 

obtained in this limit, although constant, is generally not convergent for any 
finite value of e. In fact, it is only a formal or asymptotic series with 
the property that beyond some level of approximation the series begins to 
diverge without bound (refs. 10-12). This divergence is due to small divisors 
that appear in the expansion for NpU, and occur as a result of the fact that 
Np is an unbounded operator on the space of almost periodic functions. The 
asymptotic series derived in this manner is essentially the one obtained by 
using the averaging and perturbation techniques developed in references 3 
and 4. According to the above arguments, one may show that if the first N 
terms of a convergent series expansion for B are constant (with a)p f 0) , 
then these terms are equal to the corresponding first N terms in the asymp- 
totic expansion of B (with Up = 0) . 

If the procedure described above for obtaining the matrix B(t) = const 
up to some level of approximation is inapplicable or impractical (e.g., if 
A (t) contains dominant low-frequency terms), it may still be possible to put 
B in a form that makes equation (14) easier to analyze than the original sys- 
tem. An important example occurs when B(t) or any of its related approximate 
systems have very slowly varying coefficients. In this case, it is assumed 
that A(t) has certain properties that permit B (or any approximation) to be 
expressed in the typical form 
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B = eB(e h t, e) 

where h > 2 is an integer. For systems of this type, singular perturbation 
methods, applicable to systems containing a large parameter, are available 
(refs. 1 and 12). 

In those cases where it is not possible or practical to choose the cutoff 
frequency in such a way that B is constant or a very slowly varying matrix 
to an acceptable order of approximation, then it is generally not possible to 
find a theoretical solution of the reduced system, and numerical or analog inte- 
gration methods must be used. If the original system matrix contains only low- 
frequency elements, then there is usually no advantage in using the averaging 
method. On the other hand, if A(t) contains both high- and low-frequency ele- 
ments, the numerical analysis may be greatly facilitated (i.e., speeded up) by 
first filtering out the high-frequency elements through a suitable averaging 
transformation . 

If, for the case described above, a complete solution is not sought but 
only information on stability is required, then the second method of Liapunov 
(ref. 2) can be used to advantage. In its most simple form, this theory 
states that the system y = B(t)y will be asymptotically stable if the eigen- 
values of the symmetric matrix S(t) =(1/2)(B' + B) (where B' is the transpose 
of B) are negative for all t. One should note that this is only a sufficient 
condition for stability and that the system may be stable even if the condi- 
tion is violated. 


Relationships Between Exact and Approximate Solutions 

Since it is generally not possible to obtain the exact solution for U 
(and B and V) J it is important to know the relationships that exist between 
the exact solutions of the reduced system and the approximate ones obtained by 
terminating the successive approximations or series expansions for B at some 
level. We note that if A and U (and B) are expressed in the series form of 
equations (21) then the error involved in terminating the series at some level, 
say, e^, is of order eN+1. One can also show that a similar statement is 
true for the process of successive approximations. It is thus clear that 
under appropriate conditions the behavior and stability of the exact and 
approximate solutions will be essentially the same. These conditions depend 
on the properties of the matrix A(t), the level of the approximation N, and 
the types of behavior and stability under investigation. Because of the com- 
plicated nature of these conditions, it is impossible to give a detailed 
account of them here. Only a brief and qualitative discussion will be given. 

The most important principle in relating the behavior and stability of 
the exact and approximate systems may be stated as follows: If some or all of 

the solutions to an approximate system can be shown to grow or decay exponen- 
tially (over an appropriate time interval) at an absolute rate greater than a 
quantity proportional to the error or remainder terms omitted in the approx- 
imation (e.g., £^ + l), then there will be corresponding exact solutions with 
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the same behavior. If the Nth approximation to B, say, B N , is a constant 
matrix, then the rates of growth or decay of the approximate solutions for Y 
are given by the real parts of the eigenvalues of Bjq. If the approximation, 
B^, is a very slowly varying function, for example, of the form eB^Ce^t, e), 
h > 2, then the rates of growth or decay of the approximate solutions are 
given, to order e* 1 , again by the real parts of the eigenvalues of B^. In 
this case, however, solutions may grow over some intervals and decay over 
others . 

In the event some solutions of an approximate system neither grow nor 
decay (which may occur when the approximation MpU ~ B^ is constant with some 
zero or imaginary eigenvalues), the stability of the exact solutions cannot be 
inferred except in special cases. All that can be said generally is that the 
exact solutions cannot grow or decay at a rate that is (usually) greater than 
a quantity proportional to the remainder terms in the approximation for B. 

For a more complete discussion of these questions, including special cases, 
see references 2, 5, and 13. 


Example 

In applying the results discussed in previous paragraphs, we consider the 
following second-order system 


ii + 2w 0 CU + co Q 2 [l + 2c(>(t)]u = 0 (23) 

which describes the motion of a damped harmonic oscillator parametrically 
excited by a forcing function <j>(t). It is assumed that co 0 and z, are con- 
stant and c. and <j>(t) ate small. In addition, it is assumed that <J>(t) can be 
represented as a generalized Fourier series of the form 


<P (t) 


= £<t> n cos (U) n t - e n ) 

I0 n 


By means of the transformation 

u = X-, cos co t + x 0 cos co t 
I o ^ o 

u = co 0 [-x-^ sin co Q t + x cos co Q t] 

equation (23) can be put into the form of equation (1) : 

x = A(t)x 





(24) 


(25) 


(26) 
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I 



where 


— — 


— 

— 

X 1 

A = 

a l 1 

a !2 

*2_ 


a 21 

a 22 








a il = sin 2 “o t ” ^ ” cos 2a) 0 t ^ 

a l2 = “o^ 1 " cos 2< V0 + ? sin 2a) 0 t ^ 
a 21 = + cos 2w o t ^ “ ? sin 2w o t ^ 

a 22 = -tJ o^ si - n 2u) o t + ?C1 + COS 2to Q t) ] 


( 27 ) 


To analyze equation (26) , only those properties of the solution that can be 
determined by a study of the first approximation to U, that is, U = A, will 
be investigated. The approximate system for y is then 

7 = (M p A)y (28) 

We choose the cutoff frequency of the operator M p to be some fraction 
of 2to , that is, to = 2£to 0 , l < 1. Now solutions to equation (19) for U 
exist if equation ^20) is satisfied. This is equivalent to the following con- 
dition on <(> and ?: 

|tj>| + |c| < j (3 - 2/2 )£« 0.1144 


By means of equations (11) and (27), we may write the matrix 
equation (28) as 


M A 
P 


of 


M A = 
P 


-to 


o 



+ a. 


-a 0 + a 2 

C + ai 


(29) 


where 

a 0 = Mp<f> , a x = Mp (4> sin 2u Q t) 


a 2 = M p (<f> cos 2to Q t) 


(30) 
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Using equation (24), the explicit representation of these functions is 


a o 


Z] oos^t - 0 n ) 


W <W 

n p 


"v 


a i(t) = - i Y'. 


“ n - ?y ol <w p 


4> n sin [ (w n - 2tu 0 )t - 0 n ] 


> (3D 


a 2 (t) 


1 

2 


L. 


I “n- 2w 0 1 <u p 


<j) n cos[(to n - 2co 0 )t - 0 n ] 


Unfortunately, it is not possible to solve equation (28) for completely 
arbitrary perturbation functions <f>(t) . For this purpose, (j)(t) must be fur- 
ther restricted. We will consider two cases in which a Q , a x , and a 2 are (1) 
constants and (2) very slowly varying functions. 


Case 1: Constant coefficients .- Here we assume <Kt) has no frequency 
components in the bandwidths 0 < a) n < Wp = 2 £cj 0 and | ca n - 2a> 0 1 < a)p = 2&a) 0 

with the exception of one at u> n = 0 and one at = w r = 2w 0 . In this case, 
choosing 0 O = 0, we have 

a o = ^o * a i = j <$> r sin 0 r , a 2 = \ 4>r cos e r (32) 

The frequency m r = 2a> 0 is customarily called the principal resonance 
frequency. 


The behavior and stability of the (approximate) solutions for 
determined from the eigenvalues of MpA. These eigenvalues are 


y can be 


x i = 111 ^ a i 2 + a 2 2 " a o 2 ^ 


= 0 ) 


o 



i = 1, 2 


(33) 


The approximate solutions for y(t) behave like exp^t). 

Applying the arguments made in the section concerning the relationship 
between exact and approximate solutions, one can state that the exact solu- 
tions for y will be asymptotically stable (or decay exponentially with 
increasing time) if 
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where y is a small, positive, second-order quantity (proportional to the 
squares £ 2 , <f> n 2 , etc.) and Re denotes the real part of the indicated expres- 
sion. The condition that the solutions be asymptotically unstable (or that 
they grow exponentially) may be obtained from equation (34) by reversing the 
inequality and sign of y. 

When £ = 0, a> 0 = 1, and 4> (t) is the periodic function given by 

24 > (t) = 6 - 1 + e cos 2t 

= 2(<}> 0 + <f, r cos 2t) (35) 

equation (23) reduces to the Mathieu equation. For sufficiently small values 
of <j> 0 and <f> r , or 5-1 and e, the matrix MpU of the reduced system, equa- 
tion (28), is constant. The eigenvalues of the first approximation system are 
then given by 


*i = 



(6 - l ) 2 


(36) 


The curves in the 5 - e plane separating zones of stability and 
instability are the transition curves. To a first approximation, these may be 
obtained by setting X^ = 0 in the above expression. This gives 

« = 1 ± f (37) 


which has been obtained previously by other methods (refs. 4, 5, 10, 14, 
and 15). 

Case 2: Very slowly varying coefficients.- As remarked previously, 

solutions of equation (28) cannot be found if <j> (t ) or the functions a 0 , ai , 
and a 2 of equations (31) are completely arbitrary. An exception occurs when 
these functions are very slowly varying. To be specific, we assume that all 
the frequencies appearing in these functions are small and of second order, 
that is, of order e 2 where e is of the same order as £ and <(>. In this 
case, equation (28) can be analyzed by singular perturbation methods (refs. 1 
and 12) . The approximate solutions for y have the form 

y ~ exp|J'[X i (t) + 0(e 2 )]dt} (38) 
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where A^(t) are the (time dependent) eigenvalues of the matrix MpA and are 

given, as before, by equation (33). From this it can be seen that, depending 
on whether the real parts of each A^(t) are positive or negative, solutions 
may grow over certain intervals and decay over others. 

It is of interest to estimate the conditions on the coefficients cj> 
of equation (31), for which solutions will decay monotonically (i.e., will be 
asymptotically stable) . This is the same as requiring that the real parts of 
each A^(t) be negative for all t. A sufficient condition to insure this 
type of behavior may be obtained from equations (31) and (33) and is given by 


! 

-5 ± jZXI < ' ^ (39) 

where y is a small second-order quantity and where the summation is over all 
coefficients in the series for ax (t) or a 2 (t) (i.e., eqs . (31)). 

It is interesting to compare the above stability criterion with one 
obtained by using the second method of Liapunov (ref. 2) . In this method, no 
restrictions on the coefficents a Q , ax, and a 2 in equation (28) need be made. 
In its simplest form, this theory states that the solutions to equation (28) 
will be asymptotically stable if the eigenvalues of the symmetric matrix 
S = (1/2) (B 1 + B) are negative for all t. These eigenvalues are 

A = oj 0 ±/ a i 2 + a 2 2 ] 

Taking into account the explicit form of the functions ax (t) and a 2 (t) given 
by equations (31) , it is easily shown that a sufficient condition for these 
eigenvalues to be negative is simply 
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This is essentially the same result as that given by equation (39) , which was 
based on the assumption that ax (t) and a 2 (t) are very slowly varying 
functions . 


CONCLUDING REMARKS 


In the preceding paragraphs, we have described a method of averaging 
applicable to linear, homogeneous, almost periodic differential equations, 
which is a generalization of the well-known methods developed in references 3 
and 4. Using the procedure, we have shown that it is always possible to trans- 
form an almost periodic system with rapidly varying coefficients into one with 
more slowly varying coefficients (i.e., with only low-frequency components). 

We have indicated some of the conditions under which the reduced system may be 
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more easily analyzed than the original system and have illustrated these by a 
simple example. We have also indicated a theoretical and practical advantage 
of the new method over the previous methods. This advantage stems from the 
fact that the processes used in constructing the transformation and the 
reduced system are convergent. Thus, estimates of the error incurred in term- 
inating the approximations at some level are obtained. In the previous 
methods, the processes are nonconvergent and error estimates cannot be 
obtained. Finally, we have shown how the error estimates are useful in rela- 
ting the stability and behavior of the exact solutions to the approximate 
solutions . 
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